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Abstract 



Aims. We aim at modeling the infrared galaxy evolution in a way as simple as possible and reproduce statistical properties among 
which the number counts between 15 yum and 1.1 mm, the luminosity functions, and the redshift distributions. We then aim at using 
this model to interpret the recent observations (Spitzer, Akari, BLAST, LABOCA, AzTEC, SPT and Herschel), and make predictions 
for Planck and future experiments like CCAT or SPICA. 

Methods. This model uses an evolution in density and luminosity of the luminosity function parametrized by broken power-laws 
with two breaks at redshift ~0.9 and 2, and contains the two populations of the Lagache et al. (2004) model: normal and starburst 
galaxies. We also take into account the effect of the strong lensing of high-redshift sub-millimeter galaxies. This effect is significant 
in the sub-mm and mm range near 50 mJy. It has 13 free parameters and 8 additional calibration parameters. We fit the parameters to 
the IRAS, Spitzer, Herschel and AzTEC measurements with a Monte-Carlo Markov chain. 

Results. The model ajusted on deep counts at key wavelengths reproduces the counts from the mid-infrared to the millimeter 
wavelengths, as well as the mid-infrared luminosity functions. We discuss the contribution to the cosmic infrared background (CIB) 
and to the infrared luminosity density of the different populations. We also estimate the effect of the lensing on the number counts, 
and discuss the recent discovery by the South Pole Telescope (SPT) of a very bright population lying at high-redshift. We predict 
the contribution of the lensed sources to the Planck number counts, the confusion level for future missions using a P(D) formalism, 
and the Universe opacity to TeV photons due to the CIB. Material of the model (software, tables and predictions) is available at 
http://www.ias.u-psud.fr/irgalaxies/. 
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1. Introduction 



Galaxies: statistics - Galaxies: evolution - Galaxies: star formation - Infrared: galaxies • 



The extragalactic background light (EBL) is the relic emission 
due to galaxy formation and accretion processes since the 
recombination. The infrared (8 pmi < A < 1000 fxm) part of this 
emission called cos mic infrared backgr ound (CIB) was detected 
for the first time by Pu getet alj (1 19961) and contains about half 
of the energy of the EBL (Dole et al] 120061: [Bethermin et al. 
I2010ah . Nevertheless, in the local universe, the optical/UV 
emissions are 3 times larger than infrared/sub-m illimeter ones 
(ISoife r & Neugeba uerlll99ll: iDriver et~ai1 120081 ). This pseudo- 
paradox is explained by a strong evolution of the properties of 
the infrared galaxies. 

The infrared luminosity density is dominated by 
normal galaxie s {L m . bo i ometric < 1O U L ) in the lo- 
cal Universe dSaunders et al.1 Q990). At higher redshift, 
it is dominated by luminous infrared g alaxies (L IRG, 
lO'^ o < L IR j,oiometric < 1O 12 L ) at z=l dLe Floc'h et all 
120051) and by ultra-luminous infrared galaxies (ULIRG , 
1O 12 L < L 1RMIomelric < 1O 13 L ) at z=2 (ICanuti et al.l l 2007h . 
The infrared lumin osity of these galaxies is correlated to the 
star formation rate (Ken nicuttlll998l) . Thus modeling this rapid 
evolution of the infrared galaxies is very important to understand 
the history of the star formation. 



The phy si cal m o dels (such as Lacey et al.l (HblO); 
IWilman et all (l2010h : lYounger & Hopkins! (l2oToT for the 
latest) use a physical approach based on semi-analytical recipes 
and dark matter numerical simulations. They use a limited set of 
physical parameters, but they nowadays poorly reproduce some 
basic o bservational cons traints like the infrared galaxy number 
counts (lOliveretal .1120101) . 



Th e backwards evolu t ion m o dels (like Lagache et al. 
(120041): iFranceschini et all (l2010h : iRowan-Robiiiso n (2009); 
Valian te et al . (2009)) use an evolution the luminosity function 
(LF) of the galaxies to reproduce empirically the galaxy counts, 
and other constraints. These models make only a description of 
the evolution and contain little physics. The parameters of these 
models were tuned m anually to fit observational constraints. 
iLe Borgne etail (|009) used an other approach and performed 
a non-parametric inversion of the counts to determine the 
LF. Nevertheless, this approach is complex, uses only one 
population of galaxy, and does not manage to reproduce the 
160 //m numb er counts. An other fully-empirical approach was 
used by Domingu ez et al.l (l2010l) . They fitted the SED from 
UV to mid-infrared of detected galaxies and extrapolated the 
far-infrared spectral energy distribution of these galaxies and 
the contribution of faint populations. Nevertheless, their model 
aims only to reproduce the CIB; however its ability to reproduce 
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other constraints like the number counts was not tested. 

The Balloon-borne Large-Aperture Submillimeter Telescope 
(BLAST) experiment dPascale et all l2008t iDevlin et all [2009) 
and the sp ectral and ph otometric imaging receiver (SPIRE) 
instrument |Gr iffin 2010) onboard the Herchel space telescope 



(Pilbr att & al. 1201 Oh performed recently new observations 
in the sub-mm at 250, 350 and 500 fim. In their current 
version, most of the models fail to reproduce the number 
counts measured at thes e wavelengths dPatanchon et"ai1 T2 009: 
Beth ermin etail 1201 Obi: IClements et all 1201 Ot lOliveretall 
20101) . The IValiante etal] d2009l) model gives the best results, 
using a Monte Carlo approach (sources are randomly taken in 
libraries) to simulate the temperature scatter and the hetero- 
geneity of the populations of active galactic nucleus (AGN), but 
this model strongly disagrees with the recent measureme nts of 
the redshift distribution of the CIB bv lJauzac et al.l (12010). It is 
thus necessary to develop new models that reproduce the recent 
far-infrared and sub-mm observations. 

Th e discovery o f very bright and high-redshift dusty galax- 
ies by IVieira etal] d2009h with the south pole telescope (SPT) 
suggests that the contribution of high-redshift galaxies strongly 
lensed by dark matter halos of massive low-redshift galaxies 
on the bright sub-millimeter and millimeter counts is non 
neglig ible. This contribution was discussed by iNegrello et al.l 
(2007p and an observati onal evidence of thi s phenomenon was 
found very recently bv lNegrello et alJ d2010h . We can also cite 
the simplified approach of iLima et al.l ([2010) who reproduce the 
AzTEC and SPT counts using a single population of galaxies 
with a Schechter LF at a single redshift an d a lensing model. 
We ca n also cite the very recent work of iHezaveh & Holder! 
(2010) on the effect of the lensing on the SPT counts, based on 
an advanced lensing model. 

We presen t a new simple and parametric model based on 
lLagache et al.l d2004l) SEE) libraries, which reproduces the 
new observational constraints. The parameters of this model 
(13 free parameters and 8 calibration parameters) were fitted 
from a large set of recent observations using a Monte-Carlo 
Markov chain (MCMC) method, allowing to study degeneracies 
between the parameters. This model also includes the effects 
of the strong lensing on the observations. We make predictions 
on the confusion limit for future missions, on the high-energy 
opacity of the Universe and on the effects of the strong lensing 
on the counts. This model is plugged to a halo model to study 
the sp atial distribution o f the infrared galaxies in a companion 
paper dPenin et al.l 12010). Note that an other study using also 
MCMC methods was performed bv lMarsden et al.l d2010l) at the 
same time than ours. 

W e use the WMAP 7 year best-fit ACDM cosmology in this 
paper dLarson et al.ll2010l) . We thus have H = 71 km.s^'.Mpc 1 , 
Q A = 0.734 and fi m = 0.266. 



2. Approach 

The backward evolution models are not built on physical 
parameters. Each model uses different evolving populations to 
repro duce the observational c o nstraints. Some re c ent m odels 
(like iFranceschini et al.l (120101) : iRowan-Robinsonl d2009l) ) use 
4 galaxy po pulations evolving separately to reproduce the 



on a very large library of templates and claim that the contri- 
bution of the AGNs and the dispersion of the dust temperature 
of the galaxies must be taken into account to reproduce the 
observational constraint. Our approach is to keep the model as 
simple as possible, but to use advanced methods to constrain its 
free parameters. This new parametric model can be used as an 
input for halo modeling or P(D) analysis for instance. 

As it will be shown, we did not need AGN contribution and 
temperature dispersion to reproduce the current observational 
constraints. In fact, in the l ocal Universe, the AGNs only 
domin ate the ULIRG regime (Imanishi 2009). Alexander et al] 
(2005) estimate an AGN contribu tion of 8% f or the submil- 
limeter galaxies (SMG). Recently iFadda etal] d2010l) showed 
that the proportion of AGN-dominated sources is rather small 
for LIRGs at z~ l (5%) and ULIRGs around z~2 (12%). 
Jauza c et al. showed that AGN contribution to the CIB is 

less than 10% at z<1.5. These category of luminosity dominates 
the infrared output at their redshift. The low contribution 
of AGN in these categories explains why the AGNs are not 
necessary to reproduce the mean statistical properties of the 
galaxies. Nevertheless, despite their small contribution to the 
infrared output, the AGNs play a central role in the physics of 
galaxies. 

Our model takes into account the the strong-lensing of high 
redshift galaxies by dark matter halos of elliptical galaxies. 
According to the results of Sect. 17.31 the effect of the lensing 
on the counts we fitted is smaller than 10%. The model of 
lensing does not have free parameters. It is based on WMAP- 
7-years-best-fit cosmology and on some parameters taken at 
values given by the litterature. The lensing is thus not useful to 
reproduce the current observations, but is necessary to make 
predictions at bright fluxes (>100 mJy) in the sub-mm and mm 
range, where the effects of the lensing are large. 



3. Description of the model 

3.1. Basic formulas 

The flux density S v at a frequency v of a source lying at a redshift 
z is dijogg 1999]) is 



S v = 



(1+Z)£(1 



+z)v 



AnD 2 L {z) 



(1) 



where z is the redshift, Di is the luminosity distance of the 
source, and L(\ +Z ) V is the luminosity at a frequency (1 + z)v. The 
comoving volume corresponding to a redshift slice between z 
and z+dz and a unit solid angle is 

dz V"A + (i+z) 3 a„ 

where Dh is the Hubble distance (D# = c/Hq), Da the angular 
distance to the redshift z. Q„, and £2a are the normalized energy 
density of the matter and of the cosmological constant. 

3.2. Bolometric luminosity function and its evolution 

We assume that the lu minosity function (LF ) is a classical double 
exponential function (ISaunders et al.| [T990) 



4 galaxy populations evol ving separately to reproduce me a ]r \ l ir 

observations. IValiante et al] d2009l) take randomly galaxy SEDs = X (— ) " X exp[-—log 1Q (l + — ))J (3) 
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Figure 1. Solid line: Local infrared bolometric luminosity func- 
tion from our best-fit model. Red dotted line: contribution of the 
normal galaxies. Blue dashed line: contribution of the starburst 
galaxies. Black vertical long dashed line: luminosity of the tran- 
sition between the two population (L pop ). 

where <1>(L/r) is the number of sources per logarithm of lumi- 
nosity and per comoving volume unit for an infrared bolometric 
luminosity Li R . is the normalization constant characterizing 
the density of sources. L* is the characteristic luminosity at 
the break. 1 — a and 1 - a - 1 / a 2 / In 2 (10) are the slope of the 
asymptotic power-law behavior at respectively low and high 
luminosity. 

We assume a continuous evolution in luminosity and in den- 
sity of the luminosity function with the redshift of the form 
L* oc (1 +zY L and <E>* oc (1 +z) r * where and are coefficients 
driving the evolution in luminosity and density, respectively. It is 
impossible to reproduce the evolution of the LF with constant 
and r$. We consequently authorize their value to change at some 
specific redshifts. The position of these breaks are the same for 
both ri and r$. The position of the first redshift break is a free 
parameter and converge to the same final value for initial values 
between and 2. To avoid a divergence at high redshift, we also 
add a second break fixed at z=2. 

3.3. Spectral energy distribution (SED) of the galaxies 

We use the lLagache et all (12004 SED library. This library 
contains two populations: a starburst one and a normal one. 
This library is parametrized only by the infrared bolometric 
luminosity (L/r). There is no evolution of the SED with the 
redshift. The normal population has a spectrum typical of spiral 
galaxy. The SED of this population does not evolve with L; R . On 
the contrary, the starburst SED evolves with Lj R . The brighter 
the starburst galaxy, the hotter the dust. 

The normal galaxies are dominant at low luminosity and the 
starburst at high luminosity. We thus chose arbitrary the follow- 
ing smooth function to describe the fraction of starburst galaxies 
as a function of the bolometric luminosity Lj R : 



starburst 



1 + th[logio(LiR/L pop )/o- pop \ 



(4) 



where th is the hyperbolic tangent function. L pop is the lumi- 
nosity at which the number of normal and starburst galaxies are 
equal. cr pop characterizes the width of the transition between the 
two populations. At Lj R = L pop , the fraction of starburst is 50%. 
There are 88% of starburst at L IR = L pop x and 12% at 

L/s = L pop x lO -0 "'"''. The contribution of the different popula- 
tions to the local infrared bolometric LF is shown in Fig.Q] 



3.4. Observables 

The number counts at different wavelengths are an essential con- 
straint for our model. The source extraction biases are in general 
accurately corrected for these observables. The counts are com- 
puted with the following formula 



pop 



r 

Jo 



X! I fpop(Lm) 



dN 



— (Sr) 

dSydQ. 
dLm dV 



dz 



dL IR dV^ s ^-P°P^ dS v dzdQ. 

y r_dN_ 



(5) 
(6) 

(7) 



where dN/dS v /dQ. is the number of source per flux unit and 
per solid angle. f pop (Lj R ) is the fraction of the sources of a 
given galaxy population computed with the Eq.|4] dN/dLj R /dV 
is computed from the Eq.|3] 



dN 



dN 



= = <HLir) 

dL IR dV dlog 10 (L IR )L IR log(10)dV L IR log(lOY 



(8) 



Li R (S v ,z,pop) and dLi R /dS v were computed on a grid in S v 
and z from the cosmology and the SED templates. These grids 
do not depend on the evolution of the LF nor on the population 
mixing parameters. These grids are thus generated only once 
and saved to accelerate the computation of the counts. Note that 
with this method, it is very easy to change the SED templates 
and/or add other populations. 

Other measurements help to constraint our model. For ex- 
ample, the monochromatic luminosity O mono function at a given 
redshift is 



^ fpop(LiR(vLy))<f>(L IR (vLy)) 



dlogw(L IR ) 
d(vL v ) ' 



(9) 



We do not use the bolometric LFs, because they are biased by 
the choice of the assumed SED of the sources. 

We can also compute the redshift distribution N(z) for a se- 
lection in flux S v > S vcut with 



N(z,S cut ) 



dN 
dSydz 



dSi 



(10) 



The extragalactic background due to the galaxies at a given 
wavelength is 



. =0 Js 



B v 



dN f 00 dN 

S v — — T^rrdS v dz = I S v — — —dS v (11) 

JSy=0 



-0 dSydzdQ. Js r =o dS v dil 

and can be compared to the measurements of the CIB. 
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Parameter 






Description 




Value 


a 






Faint-end slope of the infrared bolometric LF 


1.223 ±0.044 


cr 






Parameter driving the bright-end slope of the LF 


0.406 


±0.019 








Local characteristic luminosity of the LF 


2.377 ± 0.363 


j n\ /win- 

0* (z=U) (xlO 


gal/dex/Mpc 3 ) 


Local characteristic density of the LF 


3.234 


± 0.266 








Evolution or the characteristic luminosity between and Zbreak,\ 


2.931 +0.119 


r phi- k Jz 






Evolution of the characteristic density between and Zbreak,v 


0.774 


±0.196 


Zbreak,\ 






Redshift of the first break 


0.879 


± 0.052 


'L* ,mz 






Evolution of the characteristic luminosity between Zbreak,\ Zbreak2 


4.737 


± 0.301 


' phii,,mz 






P 1 1 l~i r\n nf trip r* ri Q vci ft^n ctir* H#»nci1"\f r\\ hpt~\x7Pf»Ti 7, , , gnrl -7, . _ 
J_iVUlLlLHJll Ul L11C dld.1 clULCl lallU QCllftlLy \JL UCLWCCll inbreak 1 {.break 2 


-6.246 


± 0.458 


Zbreak,2 






Redshift of the second break 


2.000 


(fixed) 


r U,hz 






Evolution of the characteristic luminosity for z> Zbreak,i 


0.145 


± 0.460 


r plii*,hz 






Evolution of the characteristic density for z> Zbreak.i 


-0.919 


±0.651 


L pop (xl0 10 L G ) 






Luminosity of the transition between normal and starburst templates 


23.677 


± 2.704 


@ pop 






Width of the transition between normal and starburst templates 


0.572 


± 0.056 



Table 1. Parameters of our model fitted to our selection of infrared observations. The errors are derived from the MCMC analysis. 



CO 
T3 




10 100 
S [mJy] 



where c is the speed of light and cr,, the velocity dispersion in 
the halo, which depends on the cosmology, the redshift and the 
mass of the halo. 

The probability P(n min ,z s ) for a source at a redshift z s to be 
magnified by a factor greater than /i m ,„ is 



d+z s f 
4nD c ( Zs ) 



FT 

Jo Jo 



dN 



P(fJ. > tU min ,Zs) 

dV 



d(log w (M))dV 



criji > fi„ 



) — dMdz 
dz 



(15) 
(16) 



where z s is the redshift of the source, D c the comoving radial 



distance, 



dN 



is the halo mass function, and 



ilV 



is the 



d(log m (M))clV iJ """" *""~"""> dz 

comoving volume ass ociated to the red shift slice dz. We use the 
halo mass function of Reed et al. (2007). 



Figure 2. Effect of the lensing on the number counts at 850 
microns. The contribution of lensed source is multiply by 10 
to underline the effect of the lensing on the counts. Red dotted 
line: counts of non-lensed sources. Green dashed line: counts of 
lensed sources. Black solid line: total counts. 

The level of the non-correlated fluctuations (shot-noise) of 
the CIB can be easily computed from our model with the equa- 
tion: 



Jo 



, dN 
'dSydD.' 



■dS^ 



(12) 



where P$n is the level of the non-correlated fluctuations and 
S v ,cut the flux limit for the cleaning of the resolved sources. 



3.5. Effect of the strong lensing on the counts 

We use a s imple strong lensing model based on iPerrotta et alJ 
(I200lll2002h . It supposes that the dark matter halos are singular 
isothermal spheres. The cross-section cr of a halo for a magnifi- 
cation fi larger than jL/ m ,„ is 



4na 2 D 



'A,ls 



(13) 



where Daj s is the angular-diameter distance between the lens 
and the source and a is given by 



4 ^ 
a = 47T— 



(14) 



The counts derived by our model take into account the fact 
that a small fraction of the sources are gravitationally magni- 
fied. The observed number counts taking into account the lensing 
(dN/dS v /d£l) lensed are computed from initial counts dS v /dz/d£l 
with 



/ dN x 

^dSydWlensed 



(S v ) 



(17) 
(18) 



r» r-dP (z) i dN fSy \ 

Jo Jfi min dfi fidS v dzdQ.\ n ' / 

Practically, this operation is performed multiplying the vector 
containing the counts for a given redshift slice by a matrix de- 
scribing the effect of lensing. This lensing matrix has diagonal 
coefficients values around 1, and small (< 10~ 3 ) non-diagonal 
terms. These non-diagonal terms describe how magnified faint 
sources affect the counts at brighter flux. The effect on the 
monochromatic luminosity function was computed in the same 
way. We chose ^„„„ = 2 which corresponds t o the limit of 
the v alidity of the strong-lensing hypothesis dPerrotta et al.l 
|2001|) . The spatial extension of the lens ed galaxies limits th e 
maximum magnification. According to Perrotta et al. (2002), 
fi max i s in the 10-30 range . We chose to use fi max =20 in this 
paper. iNegrello et ail d2007l) used ju m „,=2 and n, nax =\5. 

Fig. [2] illustrates how number counts are affected by lens- 
ing. This figure is based on the number counts predicted by the 
model at 850 [im with a probability of magnification multiplied 
by a factor 10 to better show this effect. The green dashed line 
is contribution of the lensed sources. Due to the magnification, 
the peak of this contribution is at higher flux than for non-lensed 
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sources, and due to the small probability of lensing, the peak is 
lower than for non-lensed sources. This effect of the magnifica- 
tion on the counts become non negligible when the slope of the 
counts is very steep, like in the sub-mm and mm domain. 

4. Fitting the model parameters on the data 

Our model has several free parameters. We tried to have the min- 
imum number of parameters. We determined them by fitting the 
model to published measurements of the counts and LFs. We 
used a Monte-Carlo Markov chain (MCMC) to find the best pa- 
rameters, their uncertainties, and degeneracies. We do not fit the 
measured redshift distributions, because the cosmic variance and 
the selection effects are currently not enough accurately known. 

4.1. Data: extragalactic number counts 
4.1 .1 . Data used for the fit 

We have chosen to fit the following data: 



- Spitzer MIPS counts of IBethermin et all (l2010al) at 24, 70 
and 160 pm, 

- Herschel SPIRE [Oliver et all (1201 Oh counts at 250, 350 and 
500 pm, 

- AzTEC counts of lAustermann et all d2010l) and IScott et al.l 
(2010) at 1.1 mm. 

4.1 .2. Justification of our choice 

We fit only the differential number counts since the integral 
counts are highly correlated and the correlation matrix is rarely 
estimated. 

The number counts were measured at numerous bands 
between 15 //m and 1.1 mm. We have chosen a collection of 
points. We were guided by the reliability of the measurements 
and their error bars. 

Number counts at 15 Ami based on the infrared space obser- 
vatory (ISO) data dElbaz et al.lll999tlGruppioni et alj|2002l) and 
on the Akari data (Pearson et alJ 20101 : iHo pwood et al. 2010) 
exhibit a discrepancy by a factor of about 2, and their errors do 
not include cosmic variance. The results of these papers were 
not fitted. Nevertheless, we compared a posteriori to our results 
to check consistency in Sect.P 



We fitted the Spitzer MIPS counts of Beth ermin et alJ 
(2010a) at 24, 70 and 160 pm. These points were built from the 
data of FIDEL, COSMOS and SWIRE legacy programs. The 
errors bars take into account the cosmic varianc e. These counts 
agree with previous Spitzer measurements of Papovich et al. 
(2004j): IShupe et al.l (l2008l) : lLe Floc'h et al.l d2009h:lFraver et alJ 
(2009) and Herschel measurements of iBerta et alJ (1201 Oh (in 
which the different fields are not combined). 

At 250 ixm, 350 pm and 500 pm, we fitted Herschel SPIRE 
lOliver et al .1 (120101) counts which take into account the cosmic 
variance and the deboosting uncertai nty. These coun ts agree 
with the BLAST measurem ents of iPatanchon et al.l {2009) 
and IBethermin et al.l (2010b) and the Herschel measurements 
of IClements et alJ (120101) . We chosen lOliver et all (2010) 
counts beca use Herschel data ar e better than BLAST ones 
and because IClements et al.l (1201 Oi) counts use only Poissonian 



error b ars, which coul d be la rgely underestimated. For in- 
stance, Bethermin et al ] (I2010al) estimate that the Poissonian 
uncertainties underestimate the real sample uncertainties by 
a factor 3 for counts around 100 mJy at 160 pm in a 10 deg 2 field. 

We do not fit the 850 pm observation because of the large 
discrepancies between the submilli meter common-use r bolome- 
ter array (SCUBA) observations (ICoppin e t al. 200 6]) and the 
large APEX bolometer Camera (LABOCA) ones (Weifi et al. 
120091) . We discuss this problem in the Sect. El 

We fitted th e AzT EC measurement s at 1 . 1 mm of 
lAustermann et al.l (1201 Oh and IScott et al.l (l2010h . The area 
covered by AzTEC is small compared to Spitzer and Herschel. 
We used two independent measurements of the AzTEC counts 
to increase the weight of mm observations in our fit. 



4.2. Data: monochromatic luminosity functions 
4.2.1. Data used for the fit 

We have chosen to fit the following data: 



- IRAS local luminosity function at 60 pm of Sau nders et al.l 
(I1990I) 

- Spitzer local lum inosity function at 24 pm of 
Rodighiero et al. (200§, 

- Spitzer luminos ity function at 15 pm at z=0.6 of 

Rodighiero et al. I ll2009h . 

- Spitzer luminosity fu nction at 12 pm at z=l of 
iRodighiero et all d2009h . 

- Spitze r luminosity function at 8 pm at z=2 of ICaputi et al.l 
d2007h . 



4.2.2. Justification of our choice 

We fitted some monochromatic luminosity functions. We chose 
only wavelengths and redshifts for which no K-corrections are 
needed. These observations strongly constrain the parameters 
driving the redshift evolutions of our model. 



From the iRodighiero et al ] (120091) LFs measured with the 
Spitzer data at 24 pm, we computed 3 non K-corrected LFs at 
z=0, 0.6 and 1. We used their local LF at 24 pm. At z = 0.6 and 
1, instead of using directly their results in their redshift bins, we 
combined their 15 pm LF at z=0.6 (respectively 12 pm LF at 
z=l) in the 0.45 < z < 0.6 and 0.6 < z < 0.8 bins (respectively 
0.8 < z < 1.0 and 1.0 < z < 1.4) to obtain 15 pm LF at z=0.6 
(respectively 12 pm LF at z=l). The error on a point is the 
maximum of the combination of the statistical errors of the two 
bins and of the difference between the measures in the two bins. 
The second value is often larger due to the quick evolution of 
the LF and the cosmic variance. We fitted only the points that 
do not suffer incomplet eness to avoid possi ble biases. We also 
fitted the 8 pm at z=2 of ICaputi etaTI d2007l) . 

We also fitted the local LF at 60 pm determ ined from the 
infrar ed astronomical satellite (IRAS) data (Saunders et al. 
1 19901) to better constrain the faint-end slope of the local LF. 
Due to the strong AGN contamination at 60 pm in the ULIRG 
regime, we did not fit the points brighter than 10 L Q at 60 pm. 
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Instrument 


Calibration parameter (y^) 


Calib. uncertainty 


MIPS 24 yum 


1.00±0.03 


4% 


MIPS 70 /;m 


1 06+0 04 


1°/, 


MIPS 160 /mi 


0.96+0.03 


12% 


SPIRE 250 /mi 


0.88±0.05 


15% 


SPIRE 350 /mi 


0.97±0.07 


15% 


SPIRE 500 //m 


1.17+0.1 


15% 


AzTEC 1.1 mm 


0.98+0.09 


9% 



Table 2. Calibration parameters and l-cr marginalized errors 
from our MCMC fit compared with calibration uncertainties 
given by the instrumental teams. 



4.3. Data: CIB 

The bulk of the CIB is not resolved at SPIRE wavelengths. 
We thus used the absolute measurement of the CIB level 
in SPIRE bands as a constraint for our model. We used 
the lLagache et al.l (1 19991) measurement performed on the far- 
infrared absolute spectrophotometer (FIRAS) data: 11.7+2.9 
nW.m 2 .sr _1 at 250 /mi, 6.4+1.6 nW.m 2 .sr _1 at 350 /mi and 
2.7+0.7 nW.m 2 .sr _1 at 500 //m. We assume that the CIB is only 
due to galaxies and thus neglect a possible extragalactic diffuse 
emission. 



4.4. Calibration uncertainties 

The calibration uncertainty is responsible for correlated un- 
certainties between points measured at a given wavelength 
with the same instrument. A change in the calibration modifies 
globally the number counts and the LE Assuming the "good" 
calibration is obtained in multiplying the fluxes by a factor y, 
the "good" normalized counts are obtained with S neK = yS and 
(S 2 J od dN/dS good ) = y L5 (S 25 dN/dS). The effect on the LF in 
dex per volume unit is more simple. We just have to shift the 
luminosity in abscissa by a factor y. 

We added to our free parameters a calibration parameter 
for each fitted band (see Table |5J. We took into account the 
uncertainties on th e calibration estimated by the instrumen- 
tal team in our fit (Stansberrv et al. 2007; Gordon et al. 2007; 
lEngelbracht et al.l2007HSwinvard et al.ll2010blScott et alj2010h . 



4.5. Fitting method 

To fit our points, we assumed that the uncertainties on the mea- 
surements and on the calibrations are Gaussian and not corre- 
lated. The log-likelihood is then 



N points 



log(L(8)) = Yi 



k=\ 



(m k - m modelik {Q)Y 
2<rl 



Nbond 

z 

h=\ 



(7b ~ I) 2 

2°~calib.b 



(19) 



where L is the likelihood, 8 the parameters of the model, a 
measurement, m mo d e i.k the prediction of the model for the same 
measurement, <x„, the measurement uncertainty on it, y^ the 
calibration parameter of the band b and cr ca iib t b the calibration 
uncertainty for this band. 

We used a Monte Carlo Markov chain (MCMC) Metropolis- 
Hasti ngs algorithm (IChib & Greenberd 119951: iDunklev et al. 
120051) to fit our model. The method consists in a random walk 
in the parameter space. At each step, a random shift of the 
parameters is done using a given fixed proposal density. A 



step n is kept with a probability of 1 if L(8„) > L{9„-\) or 
with a probability L{9„) / L{0„-]) else. The distribution of the 
realization of the chain is asymptotically the same as the under- 
lying probability density. This property is thus very convenient 
to determine the confidence area for the parameters of the model. 

We used the Fisher matrix formalism to determine the pro- 
posal density of the chain from initial parameters values set man- 
ually. The associated Fisher matrix is: 



N points 



dm 



model.k 



dm, 



model.k 



l 



k=\ 



36; 



1 



^^"calib.b J 



(20) 



where 8 is a vector containing the model and calibration (yb) 
parameters. The term in brackets appears only for the diagonal 
terms corresponding to a calibration parameter. We ran a first 
short chain (10 000 steps) and computed a new proposal density 
with the covariance matrix of the results. We then ran a sec- 
ond long chain of 3 00 000 steps. The final chain satisfies the 
iDunklev et all (120051) criteria (f > 20 and r < 0.01). 

5. Results of the fit 

5.1. Quality of the fit 

Our final best-fit model have a x 2 (x 2 - ~2log(L) because all 
errors are assumed to be Gaussian) of 177 for 113 degrees of 
freedom. Our fit is thus reasonably good. The parameters found 
with the fit are given in Table Q] (the uncertainties are computed 
from the MCMC). The calibration factor are compatible with the 
calibration uncertainties given by the instrumental teams with a 
X 2 of 2.89 for 7 points (see Table [2]). The results are plotted in 
Fig. 



5.2. Comparison between the model and the observed 
counts used in the fit 



The IBethermin et all (1201 Oah points fit globally well, with some 
exceptions. Our model is lower by about 15% than two points 
around 300 //Jy at 24 /mi. These two points are built combining 
the FIDEL, COSMOS and SWIRE fields. The SWIRE fields are 
shallow fields and the counts could be affected by the Eddington 
bias. We also observe a slight under-prediction o f the bright 
(Sjo > 50 mJy) counts at 70 /mi. We also plotted the lBerta et al.l 
(2010) counts at 160 /mi measured using the photodetector 
array camera and spectr ometer (PACS) on the He rschel satellite. 
These counts agree with Bet hermin et alJ(l20 10a) and ourmodel. 



Our model fi ts glob ally well the lOliver et"al ] (120101) and 
Bethe rmin et alJ (1201 Obi) counts, excepting a slight under- 
prediction of the counts between 30 mJy a nd 100 mJy at 
500 /m i. There is a mild disagreement with the IClements et al.l 
d2010 ) counts, but their errors bars do not take into account the 
cosmic variance and are thus unde restimated. We al so plotted 
the results of the P(D) analysis of iGlenn etail (l2010h . These 
points and especially the error bars must be interpreted with 
caution (see the complete discussion in lGlenn et al. (2010)). We 
have plotted the knots of the smooth and power-law models. 
They globally agree with our model. 

Our model agrees v ery well with the Az TEC counts of 



Austerma nn et al.l d2010i) and IScott etail (12010). The contribu- 
tion of the strong lensing objects to the AzTEC counts is weak 
(<10%, see Sect.O- 
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Figure 3. (a) to (f): Differential extragalactic number counts used for the fit. (h): Monochromatic LFs at different wavelengths and 
redshifts. (a) to (h): T he fitted points a r e thick er. Black solid line: our best-fit model. Black dashe d line: 1-cr range of the model, (a) 
to (c): Red diamonds: \B6\henmn et alJ d2010a|) Spitzer legacy number counts, (c): Green trianeles: iBerta et all (120 1 0|) Hersch el/PEP 
number counts, (d) to (f): Red di amonds:\Olrver et alJ d2010h HerschelfHermes number counts. Green tria neles:\G\em\ et al. | (2 010) 
Herschel/Hermes P(D) analysis. Cle ments et al. d2010l) He r schel/ ATLAS number counts. Purple cross: Bethermin et al. (2010b) 
BLAST number counts, (g): Green triangles: Scott et afl (|2010) AzTEC number co unts in the CDFS field. Green triangles: 
lAustermann et alJ (120101) AzTEC number c ounts in the SHADES field, (h): Red plus: iRodighiero et alJ (l2009h local 24 fim LF 
(not fitted points in grey). Gre en diamonds: | Saunders et alJ d 19901) local 60 pm LF (shifted by a factor 10 on the y-axis; not fitted 
points in grey); Blue tr iangles: IRodighiero et al] d2009h 157mi LF at z=0.6 (shifted by a factor 100 on the y-axis; not fitted points in 
grey). Purple squares: Rodighier o et alJ d2009h 12 pm LF at z=l (shifted by a factor 1000 on the y-axis; not fitted points in grey). 
Cyan crosses: ICaputi et alJ d2007 ) 8 jum LF at z=2 (shifted by a factor 10 000 on the y-axis). 7 
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5.3. Comparison between the model and the observed 
monochromatic LFs 



15 microns 



Our model fits well our collection of LFs dSaunders et al 
Capu ti et al.l20 07; Rodi ghiero et al.l2009l) . excepting the bright- 
est point of ICaputi et al.1 (120071) . In Fig. [3] we arbitrary shifted 
the different LFs on the y-axis to obtain a clearer plot. Our model 
underestimates the 60 fim local LF in the ULIRG regime. It is 
expected because our model do not contain AGNs and confirms 
our choice of not fitting these points (Sect. l4~2l >. 



5.4. Comparison between the model and the observed 
counts not used in the fit 

We also compared our results with the counts at other wave- 
lengths. They are plotted on Fig. |4]and[5] The l-<x region of 
the model includes the jb uncertainty of Akari at 1 5 fim (4%, 
llshihara et al] d2010l) ), PACS at 1 10 fim (ab out 10%.lBerta et alJ 
(120 lOln and LABOCA at 850 fim (8.5%, IWeiB et alJ d2009l) ). 
The uncertainty on 77, is about t he same for LABOCA and 
SCUBA (-10%, IScott etalJ (120061) ). The uncertainties on the 
model are larger at these non-fitted wavelengths because the 
correlations between the model and the calibration parameters 
are not taken into account by the fit. 



At 15 fim, the lElbaz et al I (fl999h counts from different 
fields are not compatible bet ween them, but our co unts pass 
in the cloud of points. The iGruppioni et alJ d2002l) counts 
are significantly lower tha n our model a n d oth er works. We 
marginally agree with the iPearson eta D (1201 Oh counts. The 
Hopwood et al.l (1201 Oh counts measured with Akari in a field 
around Abell 2218 are lower than our model by about 25%. 
Nevertheless, their field is very narrow and their estimation may 
suffer from cosmic variance . Finally, we well agree with the very 
recent iTeplitz et al.l d2010b measurements performed with the 
infrared spectrograph (IRS) onboard the Spitzer space telescope. 



We compa re our c oun ts to | Hacking & Houckl (|1987|) . 
Lonsdale et al. l dl990l). iRowan-Robinson et al.l (U 990). 
Saun ders et all dl990l) . iGregorich et all d 19951) and iBertin et alJ 
d 1997b at 60 fim from IRAS data. There are disagreements 
between the different observations and some error bars may be 
underestimated, but our model globally agrees with the cloud of 
points. 

We can also compare the prediction of our model with 
lBertaetalJd2010l) counts at 1 10 fim. Our model globally agrees 
with their work. Nevertheless, our model tends to be higher than 
their measurement near 100 mJy. Observations on several larger 
fields will help to see if this effect is an artifact or not. 



At 850 fim, we very well agree w i th the P(D) analysis 
of the LABOCA data of IWeifi et al.l d2009h (s ee Fig. B- 
But, t he measurements p erformed wit h SCU BA dBorvs et alJ 
2003 : IScott et all [20061 ICoppin et all l2006h and LABOCA 
(Beel en et al . 2008) are significantly higher than our model at 6 
and 8 mJy. At low flux (<2 mJy), our model agre es very well 
with t he measurement performed in lensed reg ion dSmail et alJ 
12002b iKnudsen et aD2008l: IZemcov et al]l20ld) . 



We also compa re our model predi ctions with SPT measure- 
ments at 1.38 mm (iVieira et al.ll2009l ). At this wavelength, the 
contribution of the synchrotron emission of the local radio- 
galaxies to the counts is not negligible. Nevertheless, these 
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Figure 4. (a) to (c): Differential extragalactic number counts 
not used for the fit. Black solid line: our best-fit model. Black 
dashed line: 1-cr range of the model, (a): Red diamonds: 
Elbaz et aD (fl99 9) ISO counts. Gre en triangle: Gruppioni et al.l 
(2002) ISO counts. Bl ue squares: Pearson et all ( 2010h Akari 
counts. Purple cross: iHopwood et alJ (120 10) Akari (lensed) 
counts. Cyan plus: ITeplitz et alJ d2010h Spi tz er/IRS counts. 
: lHa< 



(b): R e d diamonds: lHacking & Houckl d 19871). Lonsdale et al.l 



dl990t). lR owan-Robinso~ etalJ d 19901) . ISaunderset al.1 d!990l) 



IGregorich et al Id 19951) andlBertin etalld 19971) IRAS counts, (c): 
Red diamonds: Bert a et a fl d2010l) Herschel/PEP counts. 



sources can be separated from dusty galaxies considering their 
spectru m. We thus compare our results with the counts of dusty 
sources. Vieira etal ] d2009l) measured counts for all the dusty 
sources and for the dusty sources without IRAS 60 fim coun- 
terpart. Our model agrees with these two measurements. Fig. [6] 
shows the counts of the non-IRAS dusty sources. The 7.2% cali- 
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850 /jm 




Figure 5. Integral number counts at 850 pm. Black solid line: 
our best-fit model. Black dashed line: l-cr range of the model. 
Grey plus: IZemcov et al.l ( 20101) combined SCUBA lensed 
counts. Blue dashed line: IWeiB et al l d2009t) LABO CA P(D) 
(Schechter model). Red diamonds: Coppin et al. (2006) SCUBA 
SHADES counts. Cyan square: dBeelen et al.l 12008) LABOCA 
counts around the J2 142 -442 3 Lya protocluster. Black plus: 
iKnudsen et al.l (l2008h combi ned SCUBA lensed counts; Green 
triang les: IScott et al.l (120061) revisited SCUBA counts. Purple 
cros s: iBorvs et al l d2003l) SCUBA HDFN counts. Yellow aster- 
isks: ISmail et al.l (120021) lensed counts. 
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Figure 6. Number counts at 1.38 mm of dusty s ources with- 
out IR AS 60 pm counterpart . Black diamonds: IVieira etaP 
(2009) south pole telescope (SPT) measurements. Black solid 
line: Total contribution of S(,o < 0.2 Jy sources. Green dot- 
dashed line: Contribution of the non-lensed sources. Red dashed 
line: Contribution of the strongly-lensed sources. Dotted lines 
l-cr contours. 



bration uncertainty of SPT is taken into account in the l-cr region 
of the model. 




5.5. Comparison with the observed redshift distributions 

In Fig. [7] we compare our model predictions with observed 
redshift distributions. At 24 pm, our model over-predicts 
by about 20% the n umbe r of sources below z=l compared 
to iLe Floc'h et all (120091) observations for the selection 
S24 > 80/Jy. Nevertheless, they exclude /^ B <20 galaxies and 
their number of sources at low redshift is thus underestimated. 



Figure 7. Redshift distribution of the S24 >80 /Jy (a), 
S 24 >300 pJy (b), S 25 o >40 mJy (c), S350 >30 mJy (d), 
S500 >20 mJy (e), Ssso >3 mJy (f), Shoo >3 mJy (g) sources. 
These measurements are not fitted. Black solid line: our best- 
fit model. Black dotted line: l-cr range of the model. Grey 
solid line: our best-fit model convolve d by a gaussian o f cr 7 = 
0.125z. Purple three dot- dashed line: |Le Borgne et al.l (2009) 
model. Green dashed line: IValiante et al.l (120091) model. R ed as- 
ter isks: ILe Floc'h et all d2009l) fa. bhlChapin et al.l <l2010h (c, <L 
eh lChapman et al.l (120051) (f ) andlChapin et all (120091) (e) mea- 
surements. Blue diamonds: Rodighier o et al.l (|2009) measure- 
ments(a, b). Cyan squares: \Desai et al.l (120081) measurements (b). 
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Our model also underpredicts the number of sources at z>3. 
But, the redshifts of the z>2 sources are only moderately 
accurate (<x- ss 0.25 for t^ B >25 at z~2). Because of the 
strong slope of the redshift distribution, a significant number 
of sources measured near z=3.5 could be sources lying around 
z=3 with overestimated redshift. If we convolve our model 
with a gaussian error with cr, = 0.125z to simulate the redshift 
unce rtainties, the mo d el and the measurements agrees (Fig. |7). 
The iLe Borgne etal] ((2009) model agrees very well with the 
measu rements, excepting at z<0.5 and z>2.5. The lValiante et al.l 
(2009) model poorly reprodu ces this observati o n. Th e same 
observables was measured by Rodighiero et alj (|2009l) . Their 
results are in agreement with lLe Floc'h et all (120091) . excepting 
at z>3, where they are higher. It could be explain by a larger o~ z 
at high redshift. 

We also compare the model with the redshift dist r ibutio n 
of 5 74 > 300iuJy sources measured bv l Le Floc'h et al.l J2009), 
iRodighiero etafl d2009t) and|DesaietalJ (|2008|). These different 
measurements exhibit disagreements below z=0.5. This differ- 
ence could be explained by the removing of the brightest optical 
sources (see previous paragraph). Our model overestimates the 
number of sources at z<0.5 by a factor 2. There is a rather 
good agreements between the models and the measurements 
between z-0. 5 and z=2.5, except a small overestimation by 
Valia nte et all d2009l) near z=2. At higher redshifts, the mea- 
surements are significantly higher than the models. It could be 
explained by two reasons: an effect of the redshift uncertainties 
and the absence of AGN in our model. 



We compare with the IChapin et al.l d2010l) redshift distribu- 
tions of the BLAST isolated sources at 250 fim, 350 fim and 
500 /urn. This selection of isolated sources does not allow to 
know the effective size of the field. We thus normalized our 
model and the measured counts to have J dN/dzdz = 1. Our 
predicted redshift distribution globally fits the measurements, 
except at low z at 250 fim and 350 fim. This difference could 
be explained by the selection of isolated sources, which could 
miss sources in struct u res at low redshift. The other models 
(iLe Borgne et all 120091: IValiante et"ai1 l2009h und erpredicts the 
number of sources at low z. IValiante et al.l (I2009h also slightly 
overpredicts the number of sources at z~1.5. 

We compared the redshift distribution of the SCUBA 
sources at 850 fim with the prediction of our model. We use 
the sele ction-corrected measu rements of Chap man et al.l {2005) 
used by Mars den et al.l (2010). All the models agrees with this 
measurement. 

We also compared the prediction of our model with the 
redshift distri bution of the sources detected at 1.1 mm by 
AzTEC dChapin et al.l 12009). A significant part of the sources 
detected at this wavelength (10 over 28) are not identified, and 
the selection is not performed in flux, but in signal-to-noise 
ratio. Consequently, the normalization of the redshift distribu- 
tion is not known. We thus use the same normalization than 
for the BLAST redshift distributions ( j dN/dzdz = 1). The be- 
havior predicted by our model agrees well with the observations. 



Recently. I Jauzac et alj d2010h has measured the contribution 
of the S24 > 80 yuJy to the CIB at 70 and 160 fim as a function 
of the redshift. Their stacking analysis allows to check the total 
far-infrared (FIR) emissions of the faint sources not resolved at 
these wavelengths. Our model agrees well with their results, ex- 
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Figure 8. Differential contribution of the S24 > 80 fiiy sources 
to the CIB as a function of the redshift at 70 fim (upper 
panel) and 160 fim (lower panel). R ed asterisks: measu rement 
by stacking in the COSMOS field dJauzac et al.ll2010t) . Black 
solid line: Our model ( 1-cr limit in bla c k dotte d line). Purple 
three dot-dashed line: Le Borgn e et al ] (12009ft model. Green 
dashed line: IValiante et al.l d2009) model. Blue dot-dashed line: 
Franceschini et al. (2010) model. 



cept near z=0.5 (see Fig. |8), where their low data points could 
come from a lar ge-scale underdensity i n the COSMOS field at 
this redshift. The Le Bor gne et al.l (l2009h model o verpredicts the 
contri bution of the 24 fim sources at z>l. The Valiant eet al.l 
(2009) model does not reproduce the trend of these measure- 
ments. I^anceschinreral] ([20T(J underestimate the contribution 
of the local sources and overestimate the contribution of z~l 
sources. 



5.6. Comparison with the measured Poisson fluctuations of 
the CIB 

Table [3] summarizes the recent measurements of the non- 
correlated fluctuations of the CIB (Psn) and the predictions 
of our model. Note that Psn depends strongly on the S VtCUt , 
the flux density at which the re solved sources are cleaned. We 
agree with the measurements of Miville-Deschenes et al ] (120021) 
at 60 fim and 10 fim, lLagache et all d2007l) at 160 fim and 
IViero et al.l (120091) at 250 fim and 350 fim. We found a value 
35% lower than Viero et all (|2009j) at 500 fim. This is consistent 
with the slight under-estimation of the counts at 500 fim by 
our model. Our model is also about 40% lower than the SPT 
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wavelength reference S„„ P S N,,„ es Psn,„wM <z mo M > 



pm mJy Jy 2 .sr 1 Jy 2 .sr 1 

60 Miville-Deschenes et al. (2002) 1000 1600+300 2089+386 0.20+0.01 

90 Matsuura et al. (20101 20 360±20 848±71 0.79±0.03 

100 Miville-Deschenes et al. (2002) 700 5800±1000 7364±1232 0.38±0.03 

160 Laeache et al. (2007) 200 9848+120 10834+3124 0.73+0.06 

250 Viero et al. (2009) 500 11700+400 11585+2079 0.81+0.08 

350 Viero et al. (2009) 400 6960+200 5048+1083 1.17+0.12 

500 Viero et al. (2009) - 2630+100 1677+484 1.59+0.21 

1363 Hall et al. (2010) 15 17+2 10+3 4.07+0.24 



Table 3. Level of the non-correlated fluctuations of the CIB at different wavelengths and comparison with the predictions of the 
model. The uncertainties on the model predictions take into account the uncertainties on jt. The mean redshift <z mm /w > of the 
contribution to the fluctuations is a prediction of the model. 



measurements at 1.36 mm (lHall et alJl2010h . It could be due to 
a lack of faint sources at high reds hift in our model. We also 
disagree with Matsuura et alJ d2010h at 90 pm within a factor of 
2. Nevertheless, they cleaned all the detected sources without 
fixed cut in flux. We took their "mean" value of 20 mJy for the 
flux cut. The high sensitivity of the measurements of the flux 
cut could thus explains this difference (for instance, a decrease 
of the flux cut by 25% implies a decrease of the fluctuations of 
19%). 

We also computed the mean redshift at which the fluctuations 
are emitted with 



< z >- 



(21) 



The results are written in Table[3] As expected, the mean redshift 
increases with the wavelength. Studying the long wavelengths is 
thus very useful to probe high redshift populations. 

5. 7. Comparison with the pixel histogram of the BLAST maps 

The quality of our counts at low fluxes in the sub-mm range can 
be tested using a P(D) a nalysis (ICondonlll974t iPatanchon et al.l 
2009: lGlennetZl l2010). Without instrumental noise, the prob- 
ability density of the signal in a pixel of the map, P(D), is given 
by: ' 



exp I J R(x)e ibJx dx - J R(x)dx\ 



P(D) = J ex P^J R(x)e' wx d.\ 

where R(x) is defined by 
C 1 dN (x\ 



e- iojD doj{22) 



(23) 



This probability distribution must be convolved by the distribu- 
tion of the instrumental noise. We also subtract the mean of this 
distribution. 

We tested our model with the deepest part of the observa- 
tions of the CDFS by the BLAST team. We kept only the pixels 
of the map with a coverage larger than 90% of the maximum 
coverage. We smoothed the signal, noise and beam map by a 
gaussian kernel with the same full width at half maximum than 
the BLAST beam. This smoothin g redu ces the effect of the in- 
strumental noise Pa tanchon et al.l (120091) . The predictions of our 
model and the BLAST pixel histograms at 250 pm, 350 pm and 
500 pm are shown in Fig. [9] The uncertainties on the model pre- 
dictions take into account the BLAST calibration uncertainties 



dTruch et alJ l2009). The model agrees rather well with the data. 
Nevertheless the measured histogram is slightly larger than the 
predictions of the model, especially at 500 pm. It is consistent 
with the slight under-estimation by our model of the counts at 
500 pm (the higher the counts, the larger the histogram). The 
clustering of the galaxies (negliged in this analysis) tends to en- 
large the histogram o f about 10% and could also contributes 
to this disagreement dTakeuchi & IshiH 20041 : Patanchon et al.l 



120091: lGlenn etalJl2010h . The IValiante et alJ d2009t) model fits 
very w ell fit the BLAST pi x el his tograms. [Le Borgne ~et al.l 
( 2009|) and iFranceschini et all d2010l) overpredicts the number 
of bright pixels at 250 pm and 350 pm (S v >50mJy). It is con- 
sistent with the fact that they oyerpre dict the counts at high flux 
dOliver et alJl20TotlGlenn etal]|2010l) . 

5.8. Degeneracies between parameters 

The Pearson correlation matrix of our model is given in Tab. 
[4] We found a very strong anticorrelation between cr and 
L*(z=0) (-0.90) and between L*(z=0) and 0*(z=O) (-0.85). 
These classical strong correlations are due to the choice of 
the parametrisation of the LF. There are also very strong 
degeneracies between the evolution in density and in luminosity 
of the LF (-0.81 between and the first break, -0.67 between 
the two breaks and -0.76 after the second break). 

There are also some slight degeneracies between the calibra- 
tion factors. The Spitzer calibration parameters are correlated 
(0.68 between 24 pm and 70 pm, 0.73 between 24 pm and 
160 pm, 0.62 between 70 pm and 160 pm). The other correlation 
implying a calibration factor are between -0.6 and 0.6. 

The marginalized probability distributions of each parameter 
and the 1, 2, and 3-cr confidence regions for each pair of param- 
eters are plotted Fig.[l0] Some distributions are not Gaussian. It 
thus justifies the use to use a MCMC algorithm. 



6. Interpretation of the results 

6.1. Evolution of the luminosity function 

Our model uses a very strong evolution of the bolometric 
infrared luminosity function to reproduce the infrared obser- 
vations. The characteristic luminosity (L*) strongly decreases 
since z=2 to now. This parameter has been divided by about 
a factor 50 from z=2 to 0. The characteristic density (0*) 
increases strongly between z=2 and z=l and slightly decreases 
between z=l and now. At z>2, the model is compatible with no 
evolution in luminosity and a slight decrease of the density when 
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Figure 10. Diagonal plots: marginalized probability distributions of each parameters deduced from the MCMC. Non-diagonal 
plots: 1-cr (red), 2-cr (green) and 3-cr (blue) regions for each pair of parameters). From left to right: a, cr, L*, (f>+, r^j z , r p hi t j z , 
Zbreak,i, r LtJnz , r phi ^ mz , r Litthz , Tpu^fa, L pop and cr pop (c.f. TableQ]). 



the redshift increases. The evolution of these two parameters are 
plotted in Fig.fTTI 



We compared our results with Capu ti et al.l d2007l) mea- 
surements per fo rmed from MIPS 24 pm observations and 
Mag nelli et alJ (l2009h measurement obtained using MIPS 
70 pm observations. These two works used a stacking analysis 
to measure the faintest points. The evolutions of L* and 0* 



only marginally agrees with these two works. Nevertheless, 
they use different fixed values of cr and a and an extrapolation 
from the monochromatic luminosity t o Lir. These choic es 
could imply some biases. We found as Capu ti et al.l d2007h a 
strong negative evolution in density between z~l and z~2. They 
found an evolution in (l+z)~ 3 - 9±L0 and we found (l+z) _62±a5 . 
Nevertheless, our value is probably biased by our non-smooth 
parametrization. This evolution is discussed in details by 
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Table 4. Pearson correlation matrix for our model. The part of the matrix concerning the calibration factors is not written to save 
space. 



ICaputiet al.l d2007). 

iReddv et ail d2008l) claimed that a ~ 1.6 at z>2. But, we 
do not need an evolution of a and <x to reproduce the obser- 
vations. Nevertheless, the infrared measurements are not suffi- 
ciently deep to constraint accurately an evolution of a. 

6.2. Evolution of the dust-obscured star formation rate 

The bolometric infrared luminosity density {pm) can be de- 
duced from the bolometric infrared LF. O ur local value of p/R 
((1.05+0.05)xl0 8 L .Mpc- 3 ) agrees with lVaccari eTaP fcolO) 
measurements (1.31^'^xlO 8 Lp.Mpc y 3 ). We also agree well 
wit h measurements at higher redshift (Rodighier o et al.l (2009) 
and lPascale et al.l (|2009) (see Fig.ll2b. pm can be converted into 
star formation rate density (SFRD) using the conversion factor 
Uxl0- w Me.yr-KL^ (iKennicutti [T998D. The SFRD derived 
from o ur model agrees rather well with the Hopk ins & Beacoml 
(2006) fit of the optical and infrared measurements. 

We also determined the contribution of the differ- 
ent ranges of luminosity (normal: L/r < 10 n L Q , LIRG: 
10 11 < L m < 1O 12 L , ULIRG:10 12 < L IR < 1O 13 L , HyLIRG: 
L/r > 10 13 L o ). Between z=0 and 0.5, the infrared luminosity 
density is dominated by normal galaxies (L/« < 1O U L ). 
Their contribution decreases slowly with redshift due to the 
evolution of the LF seen in Fig.[TT] Between z=0.5 and 1.5, the 
infrared output is dominated by the LIRG. At higher redshift, 
it is dominated by ULIRGs. The HyLIRGs never dominate and 
account for some percent at high redshift. A physical cutoff 
at very high luminosity thus would not change strongly the 
infrared density evolution. 

Following our model, the number of very bright objects 
(> lO 125 L ) is maximal around z=2 (see Fig. ITTb . These 
objects could be very massive galaxies observed during their 
formation in the most massive dark matter halos. Among other 
analysis, the study of the spatial distrib ution of the galax ies will 
help to confirm or infirm this scenario (Penin et al. 2 0101) . 

Around z=l, the number of very bright objects is lower than 
at higher redshift, but the number of LIRGs is about one order 
of magnitude larger. From z=l to now, the infrared output has 
decreased by about one order of magnitude. Our model makes 
only a description of this evolution and we need physical mod- 
els to understand why, contrary to nowadays, the star formation 



at high redshift is dominated by few very-quickly-star-forming 
galaxies, when t he associated dark matter halos grew by hierar- 
chical merging (ICole et al.ll2000t lLanzoni et al.| [2005). We also 
need an explanation of the decrease of the star formation since 
z=l. The main can didates the feedback of AGNs and starbursts 
(e.g. lBaughl (120061) ) and/or the lack of gas. 

6.3. CIB SED 

The value of the CIB at different wavelengths predicted by the 
model is given in Table [5] We found a CIB integrated value 
(over the 8-1000 jjm range) of 23.7+0 .9 nW.m~ 2 .sr -1 . It agrees 
with the 24-27.5 nW.m^.sr" 1 range of lDole et af](l2006l) . 

We compared the CIB spectrum found with our model with 
the measurements (see Fig. [T3l . Our model is always higher 
than t he lower limit given by the stacking. The Marsden et al. 
d2009l) limits are very stringent. Nevertheless, they could 
be overestimated due to the contamination due to clustering 
(lBavouzejl2008t iFernandez-Conde et al.ll2010t iBefhermin et all 
1201 Obi) . Our model is compatible with the upper limit given 
by the absorption of the TeV photons by photon-photon 
interaction with the CIB (s ee Sect. I7.21>. We g lobally agree 
with the DIRBE/WH AM dLagache et all 120001) and Akari 
(Matsuura et al.l 1201 Oh absolute measurement, excepting at 
90 urn (Akari) and 100 jum (DIRBE/WHAM) where the 
measurements are significantly higher than our model. These 
measurements need an accurate subtraction of the zodiacal light 
and of the galactic emissions and an accurate inter-calibration 
between DIRBE and FIRAS. Indeed, a bad removal o f the 
zodiacal light explains this disagreement dDole et al.ll2.006h . At 
larger wavelengt hs, we very wel l agree with the FIRAS absolute 
measurements of lLagache et al.l yOOO). 

We separated the contribution of the infrared galaxies to the 
CIB in 4 redshift slices, each slice corresponding to about a quar- 
ter of the age of the Universe (Fig. [131 . Between 8 and 30 fim, 
we can see a shaky behavior of each slice due to the PAH emis- 
sion bands. The total is smoother. The < z < 0.3 dominates 
the spectrum only near 8 yum due to the strong PAH emission at 
this rest-frame wavelength. This slice, where the infrared lumi- 
nosity density is the lowest, has a minor contribution at the other 
wavelengths. The 0.3 < z < 1 slice dominates the spectrum be- 
tween 10 and 350 fim. The sub-mm and mm wavelengths are 
dominated by the so urces lying at higher redshift (z > 2, see 
lLagache et"atl (120051) 1. It is due to redshift effects that shift the 
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Figure 9. Comparison with the BLAST pixel histogram at 
250 pm (upper panel), 350 /mi (middle panel) and 500 pm 
(lower panel). (Black histogram): histogram of the values of the 
central part of the BLAST beam-smoothed map in Jy/beam. (Red 
solid line): distribution predicted by our model using a P(D) 
analysis. Our analysis does not include the clu stering. Purple 
three dot-d ashed line: | Le Borgne et alJ (l2009h model. Green 
dashed line: IValiante et alj d2009) model. Blue dot-dashed line: 
iFranceschini et al.l d2010h model. 



peak of emission around 80 fim rest-frame to the sub-mm. The 
mean redshift of the contribution to the CIB is written in Table [5] 
and computed with 



< z >= 



r m dB r > 
Jo z W dz 



(24) 



We also separate the contribution of the different infrared lu- 
minosity classes. The normal galaxies and LIRGs dominate the 
background up to 250 pm. It is compatible with the fact that 
these populations are dominant at low redshift. At larger wave- 
lengths, the redshift effects tend to select high redshift sources; 
LIRGs and ULIRGs are responsible for about half of the CIB 
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Figure 11. Evolution of the bolometric infrared luminosity 
function with the redshift. Upper panel: bolometric LF at z=0 
(solid line), z=0.5 (dot line), z= 1 (dashed line), z=1.5 (dot-dash 
line), z=2 (3-dot-dash line), z=3 (long dashed line). Lower 
panel: Evolution of the L* (red solid line) and 0* (blue dot-dash 
line) parameter as a function of redshi ft and 1-cr con fidence re- 
gion. The measurement of L* by Caput ietaLl (120 07 ) (triangles) 
using 24 pm obervations and Magnelli et alj d2009h (diamonds) 
using 70 um observat i ons ar e plotted in red. The measuremen t 
of cp* bv lCaputi etal] d2007l) (cross) and iMagnelli et all d2009h 
(square) are in blue. 



each. The HyLIRG have only a small contribution (< 10%) in- 
cluding in the mm range (Fig.[T3] bottom). 



7. Predictions 

7.1. Confusion limit 

The confusion limit can be defined in several ways. The ra- 
dioastronomers use classically a source density criteria, where 
the confusion limit is the flux cut for which a critical density of 
source is reached. The cho ice of this critical density is not trivial. 
We follow the approach of lDole et al.l (l2003). The source density 
limit Ns dc is reached when there is a probability P to have an 
other source in a k Ofwhm radius (wher e Ofwhm is t he full width 
at half maximum of the beam profile), foole et all d2003) show 
that 



N. 



SDC 



' 7Tk 2 Q 2 

°FWHM 



(25) 
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Figure 12. Evolution of the bolometric infrared luminosity den- 
sity (black solid line) as a function of the redshift. The con- 
tribution of normal galaxies (L IR < 10 11 L Q ), LIRG(10 n < 
L IR < 1O 12 L ), ULIRG(10 12 < L, R < 10° L ), HyLIRG(L /s > 
10 13 L o ) are plotted with short-dahsed, dot-dash, three-dot- 
dash, long-dash e d line respectively. The measurements of 
iRodighiero et alj d2009l) using the MIPS 24 pm data are plot- 
ted with diamonds and Pascal e et alJ (120091) ones using a BLAST 
stacking analysis with triangles. The Hopkins & Beacom (2006) 
fit of the optical and infrared measurement is plotted with a dark 
grey area (1-cr) and a light grey area (3-cr). 
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Table 5. Surface brightness of the CIB and mean 
the contribution to the CIB at different waveleng 
by the model. 



-0201 

. redshift <z> of 
;ths as predicted 



We chose P = 0. 1 and k = 0.8 following Stole et alj d2003l) . 

This source density criterion does not take into account 
the contributions of the sources fainter than the flux limit. We 
made also an estimate of the photometric confusion noise based 
on the P(D) analysis (see Sect. I5.71 i. The P(D) distribution 
in absence of instrumental noise is not gaussian and have a 
large tail at bright flux. Thus the standard deviation is not a 
good estimator of the confusion noise. We chose to compute 
the interquartile interval of the P(D) divided by 1.349. With 
this definition, the value of the confusion noise is exactly <x in 
the Gaussian case, and we are less sensitive to the bright outliers. 

These two estimators can be computed from the counts 
predicted by our model. We assume that the sources are point- 



Figure 13. Upper panel: contribution to the CIB per redshift 
slice. Black solid line: CIB spectrum predicted by the model. Red 
short-dashed line: Contribution of the galaxies between z=0 and 
0.3. Green dot-dash line: Same thing between z=0.3 and 1. Blue 
three dot-dash line: same thing between z=l and 2. Purple long- 
dashed line: Contribution of the galaxies at redshift larger than 
2. Black a rrows: Lower limits co ming from th e number counts 
at 15 u m dHopwood et al.1 l2010h and 24 /mi (Beth ermin et al 
2010a) and the stacking analysis at 70 ^m dBethermi n et al 
2010ah . 100 um, 160 um dBerta et alJl2qi0h. 250 um, 350 jum, 
500 jum dMarsden et al ll2009h . 850 um dGreve et~aT]l2009h and 
1.1 mm dScott et all 1201 Oh and upper limits coming from ab- 
sorption o f the TeV photons of [Stecker & de Jagen d 19971) at 
20 /im and Renault et al. (2001) between 5 yum and 15 jum. Black 
diamonds: iMatsuura et al.l (12 010) absolute measurements with 
Akari. Black square: lLagache et al.l d2000|) absolute m e asure- 
ments with DIRBE/WHAM. Cyan line: lLagache et al.l (l2000h 
FIRAS measurement. 

Lower panel: Contribution of to the CIB of the normal galaxies 
(red short-dashed line), LIRGs (green dot-dash line), ULIRGs 
(blue three dot-dash line) and HyLIRG (purple long-dashed 
line) and all the galaxies (black solid line). 



like. The confusion noise found for large telescope at short 
wavelength (<8 //m for a 0.85 m-diameter telescope like Spitzer 
and <35 pm for a 3.29 m-diameter telescope like Herschel) are 
thus underestimated. For this reason, we do not estimate the 
confusion levels for beam smaller than 2 arcsec. 
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Figure 14. Upper panel: l-cr confusion noise as a function of 
the wavelength for different telescope diameters. We use the 
confusion noise given b y the P(D) method (see Sect. 17.11 ) for 
this plot. Red triangles: iFraver et ail (120091) Spifeer/MIPS con - 
fusion measurements. Black diamonds: Nguvenetal. (2010): 
Herschel/SPIRE confusion noise measurements (5-cr co „f cut). 
Lower panel: Resolved fraction of the CIB by sources brighter 
than 5-cr confusion noise (fluctuations) and the source density 
limit. 

Both panel: Red dotted line: Telescope with a diameter of 0.85 m 
like Spitzer. Black solid line: 3.29 m telescope like Herschel. 
Green dashed line: 12 m telescope like Atacama pathfinder ex- 
periment (APEX). Blue dot-dashed line: 15 m telescope like the 
CSO. Purple three dot-dashed line: 25 m like the CCAT project. 
Asterisks: Transition between the source density limitation (short 
wavelengths) and the fluctuation limitation (long wavelengths). 



The Fig. [14] (upper panel) represents the conf usion noise. 
It agre es wit h the confusio n noise measured by IFraver et alJ 
(2009)) and iNguven etafl (1201 Ol) with Spitzer /MIPS and 
Herschel/SPTRE. \Wei& et alJ d2009l) estimate that the confusion 
noise for a APEX/LABOCA map smoothed by the beam is 
0.9 mJy/beam. We find 0.6 mJy/beam with the P(D) approach. 

We also compute the resolved fractio n of the CIB by 
sources brighter than the confusion limit of iDole et alJ ([2003) 
(source density criterion) and the 5-<r con f given by the P(D). 
Fig. [H (lower panel) and Table [6] [7] H [9] and [TO] summarize 
the results. The transition in the confusion regime between 
the source density limitation (short wavelengths) and the 
fluctuation limitation (long wavelengths) happens at 100 pm 
for Spitzer, 220 pm for Herschel and 1120 pm for the CSO 
(asterisks in the lower panel of Fig. FBI . For larger antennas 
below 1 .2 mm, the confusion is mainly due to the source density. 

According to these results, at the confusion limit, Herschel 
can resolve 92%, 84%, 60%, 25.9%, 9.2% and 3.3% of the CIB 
at 70 pm, 100 pm, 160 pm, 250 pm, 350 pm and 500 pm re- 
spectively. Nevertheless, due to the blackbody emission of the 



telescope (about 60 K), very long integration times are needed 
to reach the confusion limit at short wavelengths. The confusion 
limit in PACS will be reach only in the ultra-deep region of the 
H-GOODS survey. The confusion limit will probably be never 
reached at 70 pm. A telescope with the same size as Herschel 
and a cold (5K) mirror, like SPICA, could resolve almost all 
the CIB from the mid-infrared to 100 pm. A 25 m single-dish 
sub-mm telescope like the Cornell Caltech Atacama telescope 
(CCAT) project would be able to resolve more than 80% of the 
CIB up to 500 pm. 

7.2. High energy opacity 

The CIB photons can interact with TeV photons. The cross sec- 
tion between a E y rest-frame high-energy photon and an infrared 
photon with a observer-frame wavelength Aj R i nteracting at a 
redsh ift z with an angl e (and p — co s(0)) is (lHeiflerlll954 
Uauch & Rohrh^[l976l) 



o- yy (E y ,A !R ,p,z) = H(l - ^) -B 2 ) x 

•\+B 
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l-B 



with 



i- 6 A 



e th (E y ,p,z) = 



2(m e c 2 ) 



2\2 



E y {\-p)(l+zY 
hc{\ +z) 



(26) 
(27) 

(28) 
(29) 
(30) 



i-IR 



where crj is the Thompson cross section (6.65 x 10 29 m 2 ), 
m e the mass of the electron and H the Heaviside step function 
(H(x)=l ifx>0and0else). 

The optical depth t(E 7 ,z s ) for a photon observed at en- 
ergy E y and emitted at a redsh ift z s can be easily com- 
puted (iDwek & Krennrichl 120051: lYounger & Hopkins! 120101: 
iDominguez et alJfeOlOl) with 



t(E 7 ,z 



Jo 



dz 



D, 



J-l L J5u 



V^a + (i +z) 3 n„ 



(31) 



dA, R n AlR (A, R , z)(l + z) o- yy (E y , A IR ,p, z) (32) 



where n^ m (AiR,z) is the comoving number density of photons 
emitted at a redshift greater than z between Aj R and Aj R + dAj R . 
The 5 pm cut corresponds to the limit of the validity of our 
model. The number density of photons is computed with 



nA IK (J-iR,z) 



An 

hcAn 



-CB, 



CIB 



,CMB 



(33) 



where B v qib is the CIB given by our model and B v qmb is 
the brightness of a blackbody at 2 .725K corres ponding to 
the cosmic microwave background (Fixsenl [20091) . Our pre- 
dicted opacities do not take into account the absorption by 
the cosmic optical background photons (COB, A < 5 pm). 
lYounger & Hopkinsl (1201 Ol) showed that the contribution of the 
COB to the opacity is negligible for energies larger than 5 TeV. 
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Figure 15. Fazio-Stecker relation: energy at which the opac- 
ity reach a given r as a function of redshift. This plot is done 
for t = 1 (red solid line), 2 (blue dashed line) and 5 (green 
dotted line). The data p oints are th e cutoff energy of Mkn 501 
dAharonian etal.1 1 1999h. Mkn 421 dAharonian et all 120021) and 
Lac 1ES 1959+650 dAharonian et al.ll2003l) . 



7.3. Effect of the strong tensing on the number counts 

The strongly-lensed fraction is the ratio between the counts of 
lensed sources and the total observed counts. Because the slope 
of the counts varies a lot with the flux and wavelength, this 
fraction depends on the flux and the wavelength (see Fig. [ToT i. 
The strongly lensed fraction is always smaller than 2% below 
250 pm and is thus negligible. At larger wavelengths, we predict 
a maximum of the strongly lensed fraction near 100 mJy. At 
500 pm, about 15 % of the sources brighter than 100 mJy are 
lensed. This fraction increases to 40% near 1 mm. 



Ou r results can be compared with ones of Neg rello et al.l 
(2007) model. The two model predict that the lensed fraction as 
a function of the flux is a bump around 100 mJy. But, the ampli- 
tude of this bump predicted by the two models is significantly 
different. For instance, the maximum of the len sed fraction at 
500 p m is 15% for our model and 50% for the iNegrello et al.l 
(2007 ) model. The slope between 10 and 100 mJy is steeper in 
Neg rello et al.l (2007) model than in ours and is incompatible 
with the measurements dClements et alj2010tloliver et alhoiOt 
iGlenn et"aT]|2010l) . The steeper t he slope, the la r ger the lensed 
fraction. This explains why the iNegrello et al.l d2007l) model 
predicts larger lensed fraction than ours. The probability for a 
source to be lensed increases with its redshift. Differences in 
the redshift distributions of the models could also explain some 
differences in the lensed fraction. 
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Figure 16. Fraction of strongly lensed sources (magnification 
larger than 2) as a function of the flux at 350 pm (red solid line), 
at 500 pm (green dashed line), at 850 pm (blue dot-dashed line) 
and at 1.1 mm (purple three-dot-dashed line). The dotted lines 
represent the 1-cr confidence area of our model. 



We can determine up to which redshift the opacity stays 
lower than 1. We can thus define an horizon as a function of 
the energy, called Fazio-Stecker relation. We can see in Fig. [131 
that the observed energy cutoff of lo w-redshift blazars (Mk n 501 
dAharonian et alJll999h . Mkn 421 dAharonian etall 120021) and 
BL Lac 1ES 1959+650 dAharonian et al.ll2003l) ) is compatible 
with this relation. 



The Fig. [6] shows the respective contribution of the lensed 
and non-lensed sources to the SPT counts of dus ty sources 
witho ut IRAS 60 pm counterparts at 1.38 mm dVieira et al.l 
2009). According to the model, these counts are dominated by 
strongly-lensed sources above 15 mJy. These bright sources 
are thus very good candidates of strongly-lensed sub-mm galaxy. 

We made predictions on the contribution of the strongly- 
lensed sources to the Planck number counts (see Fig. [TTT i. We 
use the iFernandez-Conde et al.l (Hp08) 5-cr limits, because they 
take into account the effect of the clustering on the confusion 
noise. This effect is non-negligible due to the large beam of 
Planck. We found that the contribution of the lensed sources to 
the Planck counts is negligible in all the bands (a maximum of 
0.47 galaxies. sr~' at 550 pm). At high redshift, Planck will prob- 
ably detect more small structures like proto-clusters, than indi- 
vidual galaxies. Planck is thus not the best survey to find lensing 
candidates. Sub-mm surveys with a sensitivity near 100 mJy are 
more efficient. For instance, the Herschel- ATLAS survey should 
found 153+26 and 411+24 lensed sources with Ssoo >50 mJy 
and S 350 >50 mJy, respectively, on 600 deg 2 . 



8. Discussion 

8.1. Comparison with other backwards evolution models 

The evolution of the infrared luminosity density predicted by 
our model can be compared with the predi ction of the recent 
backw ards evolution models. We find, like iFranceschini et al] 
d2010l) . a strong increase of p IR from z=0 to z=l, a break 
arou nd z=l, and a decrea se at larger redshift. On the contrary, 
the IValiante et all (l2009t) and ILe Borgne etaD (I2009I) models 
predict a maximum of infrared luminosity density around z=2. 



As ILe Borgne etal] d2009l) and IFranceschini et alj (2010), 
we found that LIRGs dominate infrared luminosity density 
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Figure 17. Differential number counts in the Planck bands. 
These counts takes only into account the individual star form- 
ing galaxies. Black solid line: Total contribution. Green dot- 
dashed line: Contribution of the non-lensed sources. Red dashed 
line: Contribution of the strongly-lensed sources. Dotted lines 
1-cr contours. Vertic al long-dahsed line: 5-cr lim its (confu- 
sion+instrumental) of Fernan dez-Conde et al. I J2008h for a bias 
of 1.5. 



Both models found a very strong evolution in luminosity up to 
z=2 ((1 + z) 3 4 for the Valiante et al. (2009) model; (1 + z ) 2 - 9±01 
from z=0 to 0.87+0.05 and (1 + z ) 4 - 7±0 - 3 from z= 0.87+0.05 to 
2 for our model). At lar ger redshift, our mod el is compatible 
with no evolution and the lValiante et al.l d2009t) model predicts a 
slight decrease in (1 + z) . Concerning the evolution in density, 
both m odels predicts an inc rease from z=0 to z«l (in (1 + z) 2 
for the Valian te et all d2009t) model and in (1 + z )°- 8±0 - 2 for our 
model) and a decre ase at larger redshift ((1 + z)~ 15 for the 
IValiante et all (120091) model. (1 + z )- 6 2±0 5 between z=0.87±0.5 
and z=2 and (1 + z) at z>2 for our model). These two 

models thus agree on the global shape of the evolution of the 
LF, but disagree on the values of the coefficient driving it. There 
is especially a large difference on the evolution density between 
z~ and z~2. This difference could be explained by different 
positions of th e breaks. Nevertheless, the uncertainties on the 
IValiante et al] ((2009) model are not estimated. It is thus hard to 
conclude. 



IValiante et all d2009l) and iFranceschini et al.l d2010t) used 
AGNs to reproduce the infrared observations. IValiante et all 
(120091) also used a temperature dispersion of the galaxies. Our 
model reproduce the same observations using neither AGNs nor 
temperature dispersion. This show that the AGN contribution 
and the temperature scatter cannot be accurately constraint with 
this type of modeling. 



8.2. Discriminating the models: smoking guns observations? 

Although they use different galaxy populations and evolutions, 
the backwards evolution models reproduce the number counts 
from the mid-IR to the mm domain in a reasonably good way. 
It is thus important to find new observables to discriminate 
between models. 

The comparison with the sub-mm redshift distributions of 
the bright sources is a rather simple, but very discriminant ob- 
servations. For instance, the Fig. [7] shows significant difference 
of the sub-mm redshift dist ri bution s predicted by the different 
models. The Cha pin et al.l d2010l) measurements performed 
on one small field with a cut at high flux is not sufficient to 
conclude. Herschel will help to increase the accuracy of the 
measured redshift distributions and to estimate the cosmic 
variance on them. These constraints will be crucial for the next 
generation of models. 



iJauzac et al ] (120101) showed that the redshift distribution 
of the contribution of the 24 microns sources to CIB at 70 
and 160 fim (d(vB v )/dz) is also a very discriminant constraint. 
The Fig. [18] shows the d(vB v )/dz at 350 fim. The different 
models make totally incompatible predictions in the sub-mm. 
An accurate measurement of d(vB v )/dz will be thus crucial for 
the future models. 



around z=l and that UL IRG dominate a t redshift larger than 
1.5. We also found as iLe Borgne etal] §009) that normal 
galaxies dominates only up to z~0.5. 

The IValiante et al. I d2009h and our model use a similar 
parametrization of the evolution the LF which can be compared. 



8.3. Limits of our model 

Our model is a useful tool to make a first interpretation of 
the observations from the mid-infrared to the mm domain. 
Nevertheless it is biased by some structural choice in its 
construction. 
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Le Borgne et al. (2009) 

Valiante et al. (2009) 

Franceschini et al. (2010) 




Figure 18. Differential contribution of the S24 > 80 /iJy sources 
to the CIB as a function of the redshift at 350 fim. Black 
solid line: Our model (1-cr limit in bla c k dott ed line). Purple 
three dot-d ashed line: Le Borgne et alJ d2009h model. Green 
dashed line: Valiant e et alJ (2009) model. Blue dot-dashed line: 
Franceschini et al. (2010) model. 



The choice of the parameters biases the results. For example, 
we have chosen the minimal number of parameters to reproduce 
the counts. If we would used more breaks in the evolution 
in density and in luminosity, the evolutions with the redshift 
would be smoother and the errors on the predictions would be 
different. Our errors are just the statistical errors due to the 
determination of the parameter of a given model using the data. 
It does not include the uncertainty on our hypothesis on the 
evolution (like a fixed) and the biases due to our choice of 
parameters (evolution in (1 + z)' with breaks). For instance, the 
strong decrease in density between z~0.9 and z=2 is probably 
an artifact due to our choice of parametrization. In addition, our 
model of lensing is very simple and should be updated in the 
future. Nevertheless, the contribution of the lensing in the fitted 
data is low and the bias is thus negligible. 

The backwards evolution models gives a very limited 
interpretation of the data. They are only a description of the 
evolution of the statistical properties of the infrared galaxies. 
The physical processes explaining the strong evolution of these 
objects are ignored. A more complex physical approach is 
thus necessary to deeply understand the history of the infrared 
galaxies. Nevertheless, our model is very useful to make a rapid 
interpretation of new observations and predictions for the future 
missions. 



8.4. Perspectives 

Our model fit the current data with rather simple hypotheses. 
Nevertheless, the increasing accuracy of the infrared observa- 
tions will probably help us to refine it. Lots of updates will be 
possible when we will need a more complex model. 

The a and cr parameters are fixed, but we can imagine an 
evolution of the shape of the LF with the redshift. A Fisher 
matrix analysis shows that the evolution of a at high redshift 



cannot be constraint without deeper observations in the sub-mm. 
An evolution of sigma could be constraints, but is not necessary 
to reproduce the current data. 

The evolution of the parameters is very simple in the 
current version and could be updated in using more breaks or a 
smoother functional form. 

The recent observations of Herschel will help a lot to update 
the SED used in our model, and maybe enable to determine 
their evolution with the redshift. The temperature of the big 
grain and its dispersion will be measured more accurately. 
Nevertheless, this dispersion must be modeled with a limited 
number of template to authorize MCMC approach. It is one of 
the future challenge for the evolution of our model. 

Nevertheless, each refinement increases the number of free 
parameters of the model. It is important to limit the number 
of new parameters in comparison with the number of measure- 
ments. 



9. Summary 

- Our new parametric backwards evolution model reproduces 
the number counts from 15 jum to 1.1 mm, the monochro- 
matic LF and the redshift distributions. 

- Our model predicts a strong evolution in luminosity of the 
LF up to z=2 and a strong decrease in density from z= 1 to 
z=2. We predict that the number of HyLIRG is maximum 
around z=2. 

- We find that Normal galaxies, LIRG and ULIRG dominates 
the infrared output at z=0, z=l and z=2, respectively. The 
HyLIRG accounts for a small fraction (<10%) at all red- 
shifts. 

- We reproduce the CIB spectrum and predict contributions 
per redshift and luminosity slice. We found that the mid- and 
far-infrared part of the CIB are mainly emitted by the normal 
galaxies and LIRG. The sub-mm part is mainly due to LIRG 
and ULIRG at high redshift in accordance with the sub-mm 
observations of deep fields. We estimated CIB total value of 
23.7±0.9nW.trT 2 .sr 1 . 

- We estimate the fraction of lensed sources in the sub-mm as 
a function of the flux and wavelength. This contribution is 
low (< 10%) below 500 fim, but high (up to 50%) around 
100 mJy in the mm domain. 

- We predict that the population of very bright dust y galaxies 
detec ted by SPT and without IRAS counterpart dVieira et alJ 
2009) is essentially composed of lensed sub-mm galaxies. 
We also predicts the contribution of the lensed sources to the 
Planck number counts. 

- We predict the confusion limits for future missions like 
SPICA or CCAT 

- We estimate the opacity of the Universe to TeV photons. 

- Material of the model (software, tables and predictions) is 
available at http://www.ias.u-psud.fr/irgalaxies/. 



10. Conclusion 

We showed that it is possible to reproduce the number counts 
from the mid-IR to the mm domain with a rather simple 
parametric model minimized automatically. Nevertheless, other 
automatically-tuned mod els reproduce the counts with differ- 
ent redshift distributions (Le Borg ne et alJl2009l iMarsden et alJ 
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2010). It suggests that number counts only are not enough to 
build these models. Different observables are thus crucial to dis- 
criminate the different parametrization proposed by the model 
builders. These constraints are the luminosity functions, the red- 
shift distributions, the P(D) or the fluctuations. These future 
measurements and their uncertainties have to be very robust to 
be directly fit by the next generation of models. 
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Bethermin et al.: Parametric backwards evolution model 



A 


5CT conf,P(D) 


CIB fraction" 


S conf, density 


CIB fraction 6 


/jm 


mjy 


% 


mjy 


% 


24 


5.62x10 


83.1 


7.51xl0~ 2 


72.3 


70 


3.09x10" 


51.5 


2.88x10" 


48.8 


100 


1.38x10 


36.3 


1.15X10 1 


36.1 




c Q/i v |nl 


IZ.J 


a -I a w i n 1 


1 /.Z 


250 


1.06xl0 2 


3.2 


4.41x10' 


6.9 


350 


1.13xl0 2 


0.8 


3.57x10' 


3.0 


500 


9.18x10' 


0.2 


2.24x10' 


1.4 


850 


4.12x10' 


100.0 


9.25x10° 


0.7 


1100 


2.76x10' 


100.0 


6.25x10° 


0.5 



Notes. Fraction of the CIB resolved at 5-cr con f. 

{b> Fraction of the CIB resolved at the flux limit. 

Table 6. Confusion noise and resolved fraction of the CIB at 

different wavelengths for a 0.85 m telescope(5p;fzer like). 



A 


50" C onf,P(D) 


CIB fraction" 


S conf , density 


CIB fraction" 


jim 


mly 


% 


mjy 


% 


70 


7.95X10- 2 


96.4 


1.20x10"' 


91.8 


100 


5.13x10-' 


90.8 


7.75x10-' 


83.9 


160 


5.01x10° 


67.8 


5.93x10° 


59.8 


250 


1.75x10' 


25.9 


1.28x10' 


29.6 


350 


2.30x10' 


9.2 


1.28x10' 


15.8 


500 


2.08x10' 


3.3 


9.24x10° 


8.7 


850 


1.13x10' 


1.5 


3.88x10° 


4.4 


1100 


8.40x10° 


1.2 


2.66x10° 


3.5 



Notes. {a) Fraction of the CIB resolved at 5-cr con f. 

(b) Fraction of the CIB resolved at the flux limit. 

Table 7. Confusion noise and resolved fraction of the CIB at 

different wavelengths for a 3.29 m telescope (Herschel like). 



X 


5<X conf,P(D) 


CIB fraction" 


S conf , density 


CIB fraction" 




mjy 


% 


mJy 


% 


1 AH 


Z..54X 1U 


QQ 9. 


1.U4X 1U 


QQ ^ 

yy.j 


250 


3.01x10-' 


97.6 


4.48x10-' 


92.5 


350 


1.08x10° 


88.6 


1.55x10° 


74.7 


500 


1.87x10° 


66.6 


1.86x10° 


52.4 


850 


1.55x10° 


33.8 


9.70x10-' 


29.4 


1100 


1.26x10° 


26.7 


6.89x10-' 


24.1 



Notes. <a> Fraction of the CIB resolved at 5-<r co „f. 

{b) Fraction of the CIB resolved at the flux limit. 

Table 9. Confusion noise and resolved fraction of the CIB at 

different wavelengths for a 15.00 m telescope (CSO like). 



A 


5CT con f,P(D) 


CIB fraction" 


^ con f, density 


CIB fraction" 


/im 


mjy 


% 


mJy 


% 


250 


2.81X10- 2 


99.8 


1.32X10- 2 


99.1 


350 


1.57x10-' 


98.5 


2.12x10-' 


94.2 


500 


4.31x10-' 


92.6 


6.09x10-' 


79.1 


850 


5.99x10-' 


64.6 


4.62x10-' 


49.7 


1100 


5.39x10-' 


53.1 


3.46x10-' 


41.2 



Notes. <fl) Fraction of the CIB resolved at 5-o- conf . 

{b) Fraction of the CIB resolved at the flux limit. 

Table 10. Confusion noise and resolved fraction of the CIB at 

different wavelengths for a 25.00 m telescope (CCAT like). 



A 


5(T conf,P(D) 


CIB fraction" 


S conf , density 


CIB fraction" 


/im 


mly 


% 


mjy 


% 


160 


5.86X10- 2 


99.4 


5.55xl0- 2 


98.2 


250 


7.06x10-' 


94.2 


1.11x10° 


85.6 


350 


2.08x10° 


77.9 


2.57x10° 


63.2 


500 


3.05x10° 


50.0 


2.57x10° 


41.8 


850 


2.19x10° 


23.6 


1.24x10° 


22.9 


1100 


1.74x10° 


18.4 


8.74x10-' 


18.6 



Notes. (a) Fraction of the CIB resolved at 5-cr con f. 

{b> Fraction of the CIB resolved at the flux limit. 

Table 8. Confusion noise and resolved fraction of the CIB at 

different wavelengths for a 12.00 m telescope (APEX like). 
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